% ISOM_01_model
%%%%%%%%%%%%% Parameters %%%%%%%%%%%%%%%%
% clear all
%load 2good_othertime_pvspread


T=10000;
cut=500;  % Initial condition dependence
N=1;

addpath('auxiliar')
 load seed
 rng(s);
 slength=30; %Total Length of the Sample (used+Cut)
finalcut=17; % Periods Of Market access before we use (cut)



%% 1. Simulation
S=zeros(N,T+cut,N);
S_index=zeros(T+cut,N);

inv_dist=invariant(Prob);
cum_dist=cumsum(inv_dist);
load seed
rng(s);
init0=rand(N,1);
%FInd Invariant Distribution
inva = @(r) find(r<cum_dist,1,'first'); % find the 1st index s.t. r<D(i);
% Now this are your results of the random trials
rR = arrayfun(inva,init0);

for i=1:N
rng(i,'twister');
S=markov(Prob,T+cut+1,rR(i),1:length(Prob));
S_index(:,i)=S';
end


ySIM=zeros(T+cut,N);
yNSIM=zeros(T+cut,N);
KAPPASIM=zeros(T+cut,N);
epsilons_default=zeros(T+cut,N);
epsilons_debt=zeros(T+cut,N);

rng(N+1,'twister');
epsilons_default=rand(T+cut,N);
epsilon_reentry=rand(T+cut,N);
rng(N+2,'twister');
epsilons_debt=rand(T+cut,N);
PVSIM=zeros(T+cut,N);
SSSIM=zeros(T+cut,N);
for i=1:T+cut
   for j=1:N
    ySIM(i,j)=yT(S_index(i,j));
    yNSIM(i,j)=yN(S_index(i,j));
    KAPPASIM(i,j)=KAPPAS(S_index(i,j));

   end
end


%% Core Interpolations

bpSIM=zeros(T+cut,N);
bSIM=zeros(T+cut,N);
BGSIM=zeros(T+cut,N);
BGpSIM=zeros(T+cut,N);
QSIM=zeros(T+cut,N);
DefSIM=zeros(T+cut,N);
BGp_bigSIM=zeros(T+cut,N,NBG);
qSIM=zeros(T+cut,N);
DebtLimSIM=zeros(T+cut,N);
TAUSIM=zeros(T+cut,N);
BGSIM_Defaulted=zeros(T+cut,N);

%Indices
BGiSIM=NaN(T+cut+1,N);
BGpiSIM=NaN(T+cut,N);
cSIM=zeros(T+cut,N);
priceSIM=zeros(T+cut,N);
SpreadSIM=zeros(T+cut,N);
zeroBGI_ind=find(BGgrid>-0.0001,1);
BadStandindSIM=zeros(T+cut,N);

BGiSIM(1,:)=zeroBGI_ind;
bSIM(1,:)=1;%0.5;
%%
for i=1:T+cut
    for j=1:N
    bad_stand=BadStandindSIM(i,j);
    yindex=S_index(i,j);
    BGindex=BGiSIM(i,j);
    qSIM(i,j)=q(S_index(i,j));
    bgrid_cc_temp=squeeze(bgrid(BGindex,:));
    def_temp=min(max(interp1(bgrid_cc_temp,squeeze(Def_func(BGindex,:,yindex)),bSIM(i,j),'linear','extrap'),0),1);
   if bad_stand==0 || (bad_stand==1 && theta>epsilon_reentry(i,j)) % WE ARE IN GOOD STANDING or we were in bad standing but we reenter
        if def_temp>epsilons_default(i,j) %WE ENTER A DEFAULT
            BadStandindSIM(i+1,j)=1; % TOMORROW WILL BE IN BAD STANDING
            DefSIM(i,j)=true(1);            
            BGpSIM(i,j)=0;
            BGSIM(i,j)=0;
            BGSIM_Defaulted(i,j)=BGgrid(BGindex);
            BGpiSIM(i,j)=zeroBG_ind;
            BGiSIM(i,j)=zeroBG_ind;
            BGindex=zeroBG_ind;
            bgrid_def_temp=squeeze(bgrid(BGindex,:));
            bpSIM(i,j)=interp1(bgrid_def_temp,squeeze(bp_bad(:,yindex)),bSIM(i,j),'linear','extrap');
            cSIM(i,j)=max(interp1(bgrid_def_temp,squeeze(c_big_bad(:,yindex)),bSIM(i,j),'linear','extrap'),1e-8);
            priceSIM(i,j)=(1-omega)/omega*(cSIM(i,j)/yNSIM(i,j))^(1+ita);
            QSIM(i,j)=interp1(bgrid_def_temp,squeeze(Qgov(zeroBG_ind,:,yindex,zeroBG_ind)),bSIM(i,j),'linear','extrap');
            TAUSIM(i,j)=(interp1(bgrid_def_temp,squeeze(Tau_bad(:,yindex)),bSIM(i,j),'linear','extrap'));
            SpreadSIM(i,j)=0;
            DebtLimSIM(i,j)=(KAPPASIM(i,j)*(priceSIM(i,j)*yNSIM(i,j)+ySIM(i,j)))/qSIM(i,j);
            bpSIM(i,j)=min(bpSIM(i,j), DebtLimSIM(i,j));
        else
            DefSIM(i,j)=false(1);  
            BadStandindSIM(i+1,j)=0; % TOMORROW WILL BE IN GOOD STANDING
            for jj=1:NBG
                gtemp=min(max(interp1(bgrid_cc_temp,squeeze(G_big(BGindex,:,yindex,jj)),bSIM(i,j),'linear','extrap'),0),1);
                BGp_bigSIM(i,j,jj)=gtemp+ BGp_bigSIM(i,j,max(jj-1,1)); %We build the Cumulative Distribution
                if epsilons_debt(i,j)<=BGp_bigSIM(i,j,jj) && isnan(BGpiSIM(i,j))==1 % We find the spot in which the shoc crosses the cumulative
                    BGpiSIM(i,j)=jj;
                end
            end             
            BGpSIM(i,j)=BGgrid(BGpiSIM(i,j));
            BGSIM(i,j)=BGgrid(BGindex);
            bpSIM(i,j)=interp1(bgrid_cc_temp,squeeze(bp(BGindex,:,yindex,BGpiSIM(i))),bSIM(i,j),'linear','extrap');
            cSIM(i,j)=max(interp1(bgrid_cc_temp,squeeze(c_big(BGindex,:,yindex,BGpiSIM(i))),bSIM(i,j),'linear','extrap'),1e-8);
            priceSIM(i,j)=(1-omega)/omega*(cSIM(i,j)/yNSIM(i,j))^(1+ita);
            QSIM(i,j)=interp1(bgrid_cc_temp,squeeze(Qgov(BGindex,:,yindex,BGpiSIM(i))),bSIM(i,j),'linear','extrap');
            TAUSIM(i,j)=(interp1(bgrid_cc_temp,squeeze(Tau1p(BGindex,:,yindex,BGpiSIM(i))),bSIM(i,j),'linear','extrap'));
            SpreadSIM(i,j)=(1+(delta)/QSIM(i,j)-delta)-(1+rbar);
            DebtLimSIM(i,j)=(KAPPASIM(i,j)*((priceSIM(i,j)*yNSIM(i,j))+ySIM(i,j)))/qSIM(i,j);
            bpSIM(i,j)=min(bpSIM(i,j), DebtLimSIM(i,j));
        end

   else % WE ARE IN BAD STANDING and we dont enter
          BadStandindSIM(i+1,j)=1; % TOMORROW WILL BE IN BAD STANDING
            DefSIM(i,j)=true(1);            
            BGpSIM(i,j)=0;
            BGSIM(i,j)=0;
            BGSIM_Defaulted(i,j)=BGgrid(BGindex);
            BGpiSIM(i,j)=zeroBG_ind;
            BGiSIM(i,j)=zeroBG_ind;
            BGindex=zeroBG_ind;
            bgrid_def_temp=squeeze(bgrid(BGindex,:));
            bpSIM(i,j)=interp1(bgrid_def_temp,squeeze(bp_bad(:,yindex)),bSIM(i,j),'linear','extrap');
            cSIM(i,j)=max(interp1(bgrid_def_temp,squeeze(c_big_bad(:,yindex)),bSIM(i,j),'linear','extrap'),1e-8);
            priceSIM(i,j)=(1-omega)/omega*(cSIM(i,j)/yNSIM(i,j))^(1+ita);
            QSIM(i,j)=interp1(bgrid_def_temp,squeeze(Qgov(zeroBG_ind,:,yindex,zeroBG_ind)),bSIM(i,j),'linear','extrap');
            TAUSIM(i,j)=(interp1(bgrid_def_temp,squeeze(Tau_bad(:,yindex)),bSIM(i,j),'linear','extrap'));
            SpreadSIM(i,j)=0;
            DebtLimSIM(i,j)=(KAPPASIM(i,j)*(priceSIM(i,j)*yNSIM(i,j)+ySIM(i,j)))/qSIM(i,j);
            bpSIM(i,j)=min(bpSIM(i,j), DebtLimSIM(i,j));

   end
    BGiSIM(i+1,j)=BGpiSIM(i,j);
    bSIM(i+1,j)=bpSIM(i,j);
    end
end

%% Computation

bSIM=bSIM(1:end-1,:);%

% Current Accpunt Calculation
CA_SIM=bpSIM(2:end,:)-bpSIM(1:end-1,:);  % Current Account of the Private Sector

cSIM=cSIM(cut+1:T+cut,:);
ySIM=ySIM(cut+1:T+cut,:);
yNSIM=yNSIM(cut+1:T+cut,:);
S_index=S_index(cut+1:T+cut,:);
DefSIM=DefSIM(cut+1:T+cut,:); 
BGiSIM=BGiSIM(cut+1:T+cut,:);
BGpiSIM=BGpiSIM(cut+1:T+cut,:); 
BGpSIM=BGpSIM(cut+1:T+cut,:); 
bpSIM=bpSIM(cut+1:T+cut,:);
BGSIM=BGSIM(cut+1:T+cut,:);
bSIM=bSIM(cut+1:T+cut,:);
QSIM=QSIM(cut+1:T+cut,:);
KAPPASIM=KAPPASIM(cut+1:T+cut,:);
BGp_bigSIM=BGp_bigSIM(cut+1:T+cut,:,:);
epsilons_debt=epsilons_debt(cut+1:T+cut,:);
epsilons_default=epsilons_default(cut+1:T+cut,:);
DebtLimSIM=DebtLimSIM(cut+1:T+cut,:);
qSIM=qSIM(cut+1:T+cut,:);
BGSIM_Defaulted=BGSIM_Defaulted(cut+1:T+cut,:);
CA_SIM=CA_SIM(cut:T+cut-1,:);
SpreadSIM=SpreadSIM(cut:T+cut-1,:);
BadStandindSIM=BadStandindSIM(cut:T+cut-1,:);
priceSIM=priceSIM(cut+1:T+cut,:);   
TAUSIM=TAUSIM(cut+1:T+cut,:);

VSIM=zeros(T,N);
%cSIM=zeros(T,N);
bindSIM=zeros(T,N);
TransSIM=zeros(T,N);
TightSIM=zeros(T,N);
NUSSIM=zeros(T,N);
spreadSIM=zeros(T,N);
priceSIM=zeros(T,N);
%qSIM=zeros(T,N);
for i=1:T
    for j=1:N
    bgrid_temp=squeeze(bgrid(BGiSIM(i,j),:)); 
    VSIM(i,j)=interp1(bgrid_temp,squeeze(V(BGiSIM(i,j),:,S_index(i))),bSIM(i,j),'linear','extrap');

 
   
    if DefSIM(i,j)==0
        TransSIM(i,j)=QSIM(i,j)*(BGpSIM(i,j)-(1-delta)*BGSIM(i,j))-delta*BGSIM(i,j);
        spreadSIM(i,j)=(delta-delta*QSIM(i,j))/QSIM(i,j)-rbar;
        
    else
        TransSIM(i,j)=0;
    end
    priceSIM(i,j)=(1-omega)/omega*(cSIM(i,j)/yNSIM(i,j))^(1+ita);
    DebtLimSIM(i,j)=(KAPPASIM(i,j)*(priceSIM(i,j)*yNSIM(i,j)+ySIM(i,j)))/qSIM(i,j);
   if abs(bpSIM(i,j)-(1/qSIM(i,j))*KAPPASIM(i,j)*(priceSIM(i,j)*yNSIM(i,j)+ySIM(i,j)))<1e-3
            bindSIM(i,j)=1;
            TightSIM(i,j)=bpSIM(i,j)-(1/qSIM(i,j))*KAPPASIM(i,j)*(priceSIM(i,j)*yNSIM(i,j)+ySIM(i,j));          
   end

    end
end
CSIM=(omega*cSIM.^(-ita)+(1-omega)*yNSIM.^(-ita)).^(-1/ita);
meanY=mean(priceSIM.*yNSIM+ySIM);
Y_SIM=priceSIM.*yNSIM+ySIM;
CAo_SIM=CA_SIM./Y_SIM;  % Current Account
%% Computing Solutions
 
 mean_default=mean(DefSIM);
cte_debtfactor=delta/(1-(1-delta)/R) ;
amean_PubDebt=cte_debtfactor*mean(BGpSIM)/mean(priceSIM.*yNSIM+ySIM);
amean_PvdDebt=mean(bpSIM)/mean(priceSIM.*yNSIM+ySIM);
% Private Crisis
BQSIM=(BGpSIM-(1-delta)*BGSIM).*QSIM;
avSpr=mean(spreadSIM); 

PvDebt_to_GDP=bpSIM./(priceSIM.*yNSIM+ySIM);
PubDebt_to_GDP=cte_debtfactor*BGpSIM./(priceSIM.*yNSIM+ySIM);
avPvD=mean(PvDebt_to_GDP);
avPub=mean(PubDebt_to_GDP);
 NFA=PubDebt_to_GDP+PvDebt_to_GDP;
a_Std_Spr=std(spreadSIM);

a_Std_PvDebt=std(PvDebt_to_GDP);
a_Std_PubDebt=std(PubDebt_to_GDP);
aa_pars=[kost0,kost1,meank,beta, sigmak, v];
aa_res=[avSpr,a_Std_Spr,avPvD,avPub, a_Std_PvDebt, a_Std_PubDebt];
aa_res1=[avSpr,a_Std_Spr,avPvD,avPub, a_Std_PvDebt, a_Std_PubDebt,mean(NFA),std(NFA)];

% LARGE ONES


aa_res_bis=[mean(spreadSIM),std(spreadSIM),mean(bpSIM./(priceSIM.*yNSIM+ySIM)),mean(cte_debtfactor*BGpSIM./(priceSIM.*yNSIM+ySIM)), std(bpSIM./(priceSIM.*yNSIM+ySIM)), std(cte_debtfactor*BGpSIM./(priceSIM.*yNSIM+ySIM))];

%SMall ones

aa_res_tris=[avSpr,delta/std(QSIM),mean(bpSIM)/meanY,cte_debtfactor*mean(BGpSIM)/meanY,std(bpSIM)/meanY,std(BGpSIM)*cte_debtfactor/meanY];
%% Current Account and SS
SS=zeros(T,1);
SS(CAo_SIM<(mean(CAo_SIM)-2*std(CAo_SIM)))=1; 

mean(SS)



